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ABSTRACT 

A description of the FNWC five level baroclinic, global, primitive 
equation model is presented including the finite difference equations 
used* Analytic initial data is generated in order to avoid the compli- 
cations caused by real data with respect to verification. A comparison 
of second order differencing versus mixed second and fourth order 
differencing shows that the latter gives relatively more accurate 
propagation for smaller scale disturbances than the former. In addi- 
tion various experiments are conducted to improve the model's treatment 
of cross-polar flow. 
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I. INTRODUCTION 



After the FNWC Northern Hemispheric Primitive-Equation Model became 
operational (Kesel and Winninghoff, 1972), Dr. F. J. Winninghoff, a 
former Naval Postgraduate School faculty member, designed and programmed 
a five-level global primitive equation model utilizing a spherical, 
staggered grid and the sigma coordinate system. This global atmospheric 
model was basically a conversion of the operational hemispheric model 
to a global model using a staggered latitude- longitude grid instead of 
a polar stereographic projection. 

The advantage of a global model is that it treats the entire atmos- 
phere as a physical system* It is reasonable to expect that when the 
model is fully operational, increased accuracy and range of prediction 
will result. A global atmospheric model allows tropical circulation 
to develop freely without artificial boundary conditions. This is crucial 
since the tropics are the main source region of the energy, which drives 
the general circulation of the atmosphere. Interactions between the 
hemispheres are also represented. 

The operational use of a global atmospheric model for short range 
prediction is practical since advanced computers capable of integrating 
a global model in a reasonable amount of time are presently available 
and since satellites in the near future are expected to provide more 
accurate world wide data input. 

From a tactical standpoint an operational global atmospheric model 
is a tremendous advantage since it provides a timely forecast of the 
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weather anywhere in the world. Improved predictability in the tropics 
would greatly benefit the Navy since a large percentage of its bases 
are located in the tropics and many naval operations take place in 
tropical waters. The most recent example of this is the establishment 
of U. S. Naval presence in the Indian Ocean and the building of a naval 
base on Diego Garcia Island in that area. 

Elias (1973) modified and carried out further experiments with the 
Winninghoff global model on the 6500CDC computer located at FNWC. He 
generated geopotential fields by inserting an analytic spherical harmonic 
stream function (Neamtan, 1946) into the linear balance equation and 
solving for the geopotential by successive over-relaxation. The geo- 
potential fields were generated on a 63 x 63 polar stereographic grid 
and were interpolated onto a 5° Northern Hemisphere latitude-longitude 
grid after insertion into the global model. Once this interpolation 
had taken place the fields were reflected into the Southern Hemisphere. 

Real data analyzed by FNWC objective schemes were also inserted into 
the global model on a polar stereographic grid and then interpolated 
onto a spherical grid and reflected into the Southern Hemisphere. In 
the cases where analytic data were generated, artificial moisture and 
temperature fields were also produced. The purpose of using an analytic 
geopotential field instead of real data is that it allows wave number, 
phase speed and wave amplitude to be specified. 

The geopotential fields were used in the linear balance equation to 
obtain the stream function which provided initial rotational winds. 

In the analytic cases Elias also avoided numerical balancing by using 
an analytic solution for the wind fields. His forecasts using analytically 
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derived winds were well behaved but those using winds from the linear 
balance equation excited spurious inertial-gravity waves which were 
undesirable for operational forecasts. 

Mihok (1974) and McCollough (1974) converted the FNWC global atmos- 
pheric model to operate on the IBM 360 computer at the Naval Postgradu- 
ate School. McCollough examined various schemes of initialization of 
the model using real data from the FNWC objective analyses used by 
Elias. He found that the use of the Robert (1965) time frequency fil- 
ter and dynamic balancing using forward and backward averages about 
the initial value was effective in reducing inertial-gravity wave 
"noise.” 

McCollough (1974) pointed out that the method of data input by 
interpolating to latitude-longitude grid points and to sigma surfaces 
from a polar stereographic grid on p-surfaces probably introduces some 
imbalance between mass and wind fields. In order to avoid the 
problems associated with initializing with real data this study employs 
a simplified direct analytic solution to the non-linear balance 
equation (Phillips, 1959). This solution is based on the same spheri- 
cal harmonic stream function (Haurwitz, 1940, Neamtan, 1946) used by 
Elias (1973). 

McCollough also noted that large gradients developed near the 
pole when real data were used. This thesis examines the polar problem 
with an analytic initial state which gives flow over the pole. 

Mihok (1974) modified the horizontal spatial difference equations 
to form a mixed second and fourth order scheme with the aim of improving 
the phase speeds of the meteorological waves, as proposed by Williams 
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(1972), He generated analytic fields using the same method as Elias 
(1974) and also used the same real data as Elias (1973) and McCollough 
(1974), Mihok derived his winds from the stream function computed by 
numerically solving the nonlinear balance equation with geos trophic 
winds used in the Jacobian. The use of this method of balancing intro- 
duced inertial-gravity waves which made it difficult to compare the 
relative merits of forecasts made with the second and fourth order 
scheme. The experiments in this thesis avoid the balancing problem by 
employing analytic height and wind fields. Mihok's mixed scheme is 
compared with the second order scheme for various wave numbers. 
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II* BAROCLINIC PRIMITIVE EQUATION MODEL 



The five-level global atmospheric primitive equation model currently 
under development by FNWC and the Naval Postgraduate School is used in 
this study. The primitive equations are written in spherical coordinates 
with sigma used vertically. They are similar to the equations used by 
Smagorinsky et al (1965), Arakawa et al (1969) and Kesel and Winninghoff 
(1972). Finite differencing of these equations takes place on a grid 
staggered both horizontally and vertically. 

The finite difference form of the primitive equations contain Ara- 
kawa 1 s (1966) technique for conserving energy. Fictitious kinetic 
energy is not produced by the integration of the nonlinear advective 
terms. However, difference form of the hydrostatic equation used in 
the model does not meet the restrictions placed on differencing in the 
vertical if total energy is to be conserved. A discussion of the re- 
quirements for total energy conservation may be found in Haltiner (1971). 
The difference equation used in the model is similar to the one used 
by Kesel and Winninghoff (1972). They state that the reason for 
differencing the hydrostatic equation in this manner is to reconcile 
the initial temperature, height and surface pressure analyses with the 
forecast variables in the model. 

A complete list of the finite difference equations, similar to the 
list given by Mihok (1974), is provided in Appendix A. The fourth 
order modifications to the difference equations (Mihok, 1974) are stated 
in Appendix B. 
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The heating and moisture source terms, lateral diffusion terms and 



frictional stress terms are not included in Appendix A but are the 
same as those used in the operational FNWC Northern Hemispheric Primi- 
tive-Equation Model (Kesel and Winninghoff, 1972). For the purpose of 
this study the heating and moisture source terms and lateral diffusion 
terms are turned off. 

A flat earth (all terrain heights at sea level) is used in this 
study since it was desirable to avoid the instability experienced by 
McCollough (1974) when terrain height was included. 

A. PRIMITIVE EQUATIONS 

The continuous form of the equations used in the model are as follows: 
Longitudinal momentum equation 

9(utt) _ 1 r 9 uutt 9(uvtt cos 0) 



+ + Z UV J # Xj. + •fl-yf 

9a a 



a cos 9 



1 




( 1 ) 



Latitudinal momentum equation 



9 (vtt) 1 , 9 (uvtt) 9(wrr cos 9) 

9t a cos 9 9X 99 



tt9 (wv) 
9a 



ttuu tan 9 



a 




( 2 ) 
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Thermodynamic equation 



3ttT = _ 1 

3t a cos 0 



r 3(iruT) 3 (ttvT cos 0), tt3(wT) 
L 3X 30 1 + 3a 



, RT r ^ , 1 r 3 tt . 3tt 

+ 7 T~r -WTT + a (— H [u TvT + v cos 0] ) 

C a L 3t a cos 0 3 A 30 'J 



+• TrH + D„ 



( 3 ) 



Moisture continuity equation 



3qTr 

3t 



1 r 3(TTuq) 3(TTvq) , 

a cos 0 L 3X 30 



+ 7r ^ W ^ + 1TQ + Dq 



(4) 



Mass continuity equation 



3tt 

3t 



1 r 3(u7r) 3(vtt cos 0) .. 3w 

a cos 0 1 3X 30 J 3a 



(5) 



Hydrostatic equation 



3£ = _ RT 

3a a 



(6) 



Equation of state 



aP =* RT 



(7) 



The dimensionless vertical coordinate sigma (a) is defined as 



a = — . 

TT 

The measure of the vertical motion 



( 8 ) 



3a • 
w= -_=_a 



(9) 
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A complete list of symbols and abbreviations used in the above equations 



may be found in the front of this thesis. 

The terrain pressure derivatives appearing as part of the pressure 
gradient force terms in the momentum equations are computed by a special 
finite difference scheme in order to reduce the truncation error noted 
by Kurihara (1968), Even though the main purpose of this scheme is to 
correct problems occurring over steep terrain, it was allowed to remain 
in the model during the experiments reported in this thesis in which a 
flat earth (all terrain heights at sea level) was specified. The tech- 
nique involves the expression of the gradient of terrain pressure as a 
function of geopotential on pressure and sigma surfaces. Equation (10) 
illustrates this relationship. 



The subscript p represents a pressure surface and the subscript a 
represents a a surface. The finite difference scheme given in Appendix 
A utilizes this relationship by interpolating the value of the geo- 
potential at adjacent grid points to the pressure surface of the com- 
putational grid point and then performing the integration on a pressure 
surface. The computational grid point is common to a sigma surface and 
the pressure surface on which the integration takes place. Thus, pres- 
sure force is computed at a point on the sigma surface. 

The terrain pressure derivatives appearing in the thermodynamic 
equation as part of the term representing the work done per unit time 
and mass by the pressure force are also handled in a special way as 
shown in Appendix A. 



Vrr = — [■ 



a ‘'cos 0 9X 



1 r 1 3 tt-> 




( 10 ) 
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B. SPHERICAL COORDINATES, STAGGERED GRID 

The spherical grid points occur at five degree intervals of latitude 
and longitude. These grid points are subdivided into two categories; 
mass points at which the mass variables IT, <J>, T, q, and w are retained 
and velocity points at which the motion variables u and v are carried. 
The poles are considered mass points. Mass points and velocity points 
are staggered throughout the grid, hence mass and velocity points occur 
at ten-degree intervals. Figure (1) taken from Mihok (1974) illustrates 
the spherical coordinate, staggered grid. 

When a variable is needed at a grid point on which it is not 
carried, an average of the values at the surrounding grid points is 
taken. For example if 7T is needed at a velocity point, the average 
of the 7T values at the four mass points surrounding the velocity point 
is used. 

C. VERTICAL LAYERS 

The atmosphere is vertically divided into five layers in sigma. 

The variables u, v, T, <p are retained at the 0.9, 0.7, 0.5, 0.3 and 0.1 
sigma levels, which are the centers of each of the five layers. The 
specific humidity (q) is retained only in the lower three layers. This 
is reasonable for the distribution of moisture in the atmosphere. The 
quantity w = - O is carried at the interfaces of the vertical layers 
(0.8, 0.6, 0.4, 0.2 sigma levels) and is zero at the upper and lower 
boundaries. Figure 2 provides a graphic description of the vertical 
atmospheric layering used in the model. 
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D . SPATIAL FINITE DIFFERENCE METHODS 



The experiments in this thesis fall into two categories: (1) experi- 
ments comparing second order spatial differencing versus mixed second 
and fourth order spatial differencing, and (2) experiments examining the 
reaction of the model to flow over the poles with second order differenc- 
ing. 

1. Second Order and Mixed Second and Fourth Order Finite Differencing 
Williams (1972) has shown that phase speed errors with fourth 
order finite differencing are considerably less than phase speed errors 
occurring with second order differencing. However, fourth order differ- 
encing requires a reduction in the length of the time step in order to 
insure linear computational stability. (Time step requirements are dis- 
cussed more fully in Section II-E.) Williams found that the time step 
would not have to be reduced as much if fourth order differencing is 
limited to the advection and divergence terms. These terms are the prin- 
cipal ones which affect the phase speeds of the meteorological waves. 

In order to gain improved phase speeds for meteorological waves 
while keeping the necessary reduction of the time step to a minimum, 

Mihok (1974) inserted a mixed second and fourth order differencing scheme 
into the model. The mixed scheme is documented in Appendix B. 

Another factor affecting the phase speeds of waves in the model 
is the arrangement of points in the horizontal. Phase speed errors for 
waves moving along the coordinate axes, with the difference scheme of 
this thesis, are the same as those for an unstaggered grid. However, 
Williams (1972) has noted that waves propagating at an angle of 45° to 
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the coordinate axis will have a larger phase speed error on the staggered 
grid than on the unstaggered grid. This is because the staggered grid 
points used in the difference equations are more widely separated rela- 
tive to a wave that is moving at a 45° angle with respect to the coordi- 
nate axes. 

2* Polar Finite Differencing 

One of the difficult problems associated with spherical coordinates 
is the treatment of the poles which are singular points. The zonal and 
meridional wind components are undefined at the poles. In the difference 
equations for this model any term requiring a velocity component from the 
poles is set to zero. 

Since the poles are considered mass points in this model separate 
calculations of the mass point parameters (tt, T, q, <p, w) are made for 
the poles. Following the procedure of Arakawa (1974), polar mass point 
parameters in the model change as the result of the meridional mass flux 
at all velocity points on the 85° latitude circle. The flux at each mass 
point is weighted to represent a triangular area with apexes located at 
the pole and at the mass points adjoining the velocity point. When an 
individual mass point parameter has been calculated for each of the 36 
triangular areas surrounding the pole, an average is taken and inserted 
as the pole value. As an illustration, the equation for the contribution 
of one atmospheric layer to the local change in terrain pressure can be 
written 

,9tt n 1 Autt A(v7T cos 6) nn 

~ a cos 6 " 2AA 2A0 1 ' 

where < > indicates a vertical summation of the layer contributions 
and 7T the terrain pressure at a velocity point obtained by averaging 
the tt values at the four surrounding mass points. 
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By clearing all terms from the denominator of (11) and multiply- 
ing by the earth’s radius, a , one can obtain 
♦ 

(g^*) [a 4AAA0 cos 0] =» < Autt a2A0 + A(vrr cos 0)a2AA > (12) 

2 

The quantity on the left hand side, [a 4AAA0 cos 0], represents the 
approximate area of a ten-degree square with a velocity point at its 
center. The u and v values at the central velocity point control the 
flux over the area. The triangular flux areas surrounding the pole are 
one quarter the size of the square flux areas at lower latitudes. Thus, 
the equation for each polar triangular flux area becomes 

(|^) [a 2 AXA0 cos 0] = [+ < vtt cos 0 > a 2 AX] (13) 

0 1 

All mass flux terms except the meridional flux along the 85° latitude 
circle are omitted . 

The following equation shows how the contribution of each of 
the 36 flux triangles is averaged to obtain the final value for the 
local change in terrain pressure at the pole. 

<K> ' m Te E ± <"> ' <14) 

n=*l 

An experiment was also performed using values at the poles ob- 
tained solely from averaging the desired quantity around the 85° latitude 
circle. This method is much simpler to apply and requires less computer 
time. 

Another problem which occurs near the polar region is maintaining 
computational stability as the meridians converge to the poles and the 
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grid distance along a latitude circle steadily decreases. Since it is 
not desirable to reduce the time step or modify the grid Arakawa’s (1974) 
method of smoothing the zonal component of the pressure gradient force 
and the divergence is used. A stability coefficient, given in the 
following Equation (15) is used to determine the smoothing. 



where 



S 



D cos 8 

C At sin (md) 
o 



(15) 



S = stability coefficient 
D = 5° grid distance at equator * 

C Q =* phase speed of the fastest gravity wave 
At = time step 
m = wave number 
d = 5° in radians 

If the value of S is less than 1 smoothing is necessary. The field 
to be smoothed is expanded into a Fourier series and the amplitude of 
each unstable wave component is reduced by a factor proportional to S . 
Arakawa states that this method of smoothing does not smooth the fields 
of variables because it is simply a generator of multiple point differ- 
ence quotients. 

Several of the experiments with flow over the pole involved 
modification of the stability coefficient or direct smoothing of the main 
variables in addition to the zonal pressure gradient and divergence fields. 
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E. TIME FINITE DIFFERENCE METHODS 



ward 



The initial time step in the model is performed using the Euler back- 
(Matsuno) difference scheme which is expressible as follows: 

* t 8F 1 " 

F 55 F + At [first approximation] 

* . ( 16 ) 



^t+1 t . A 8F i 

F = F + At [forecast] . 



F is a vector including the dependent variable, the superscript t 
denotes the time step, the superscript * indicates the result of the 
intermediate step and At is the time interval of a single time step. 

The Euler backward scheme is a two step iterative scheme which is either 
damping or neutral depending on the wave number (Haltiner and Williams, 
1973). The principal damping takes place in short waves. For this 
reason the Euler backward scheme is applied at three-hour intervals in 
addition to the initial application. 

All other time steps are performed using central differences, commonly 
called the leapfrog scheme; 



F 



t+1 



t-1 



+ 2At 




(17) 



The leapfrog scheme gives rise to a spurious wave called the computational 
mode which has no counterpart in the analytic solution to the differ- 
ential equation. The phase of the computational mode changes every time 
step causing the solutions at even and odd time steps to decouple. The 
amplitude of the computational mode is greatest for short wavelengths 
(Haltiner, 1971). 
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In order to limit the effects of the computational mode and reduce 
noise in the short wavelengths the Euler backward scheme is applied 
every three hours. As previously mentioned, this tends to dampen short 
waves somewhat. In addition, it is important to note that there is 
no computational mode associated with the Euler backward scheme and 
therefore application of this scheme recouples the solution. 

A time step of ten minutes is used when second order differencing 
is applied throughout the model. When the mixed second and fourth 
order differencing scheme (Williams, 1972), (Mihok, 1974) is used the 
time step is reduced to eight minutes. Williams showed, experimenting 
with the linearized barotropic primitive equations, that a 14% smaller 
time step was needed with mixed second and fourth order spatial 
differencing in order to meet the von Neumann necessary condition for 
stability. Both the ten-minute and eight-minute time steps are too 
large near the poles to maintain stability unless Arakawa f s method of 
smoothing derivatives is used. 
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III. INITIALIZATION 



Initial conditions for all experiments in this thesis were derived 
analytically. They approximate an initial barotropic state. Generation 
of initial conditions in this study differs from the method used by 
Elias (1973) and Mihok (1974) in that an analytic solution for c p rather 
than a numerical solution of the balance equation is used. Generation 
of initial conditions in this manner was found to be less time consuming 
and produced waves with almost no noise that propagated with small dis- 
tortion. 

Initial temperatures are generated on twelve constant pressure 
levels distributed from 1000 mb to 50 mb and height values are generated 
for ten of these levels. In addition, sea level pressure and sea surface 
temperature fields are produced. Wind fields are generated at seven 
constant pressure levels. All generated fields are interpolated to 
sigma surfaces. Initial values at the poles are averages taken around 
the 85° latitude circle on sigma surfaces. 

A. ANALYTIC BALANCING 

The initial geopotential and velocity fields used for the experiments 
in this thesis were obtained from the stream function solution to the 
linearized vorticity equation (Haurwitz, 1940) .Neamtan (1946) later 
solved the complete vorticity equation and obtained results verifying 
Haurwitz f s work. 
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The expression for the stream function, <p , is given by 



* 2 
if/ = A sin (mA-Vt) sin 0 cos m0 - B a sin 0 (18) 

where m is the wave number, V is the angular wave velocity, a is 

* 

the radius of the earth and A and B are constants. The amplitude 
* 

constant, A , (Elias, 1973) is defined as 

A* =* A [ ( 2m ! / 2 N N!)] , (19) 

where N =* mfl and A is a constant. The constant B is directly 
related to the angular phase speed by the equation 

V _ N(N+1) - 2 2Q 

m N(N+1) N(N+1) K } 

where (v/m) is the angular phase speed, and Q is the angular velocity 
of the earth, 

Haurwitz (1940) has shown that harmonic waves defined by this stream 
function in a barotropic non-divergent atmosphere will move with a 
constant angular velocity without changing shape. Harmonic waves 
generated in this way provide excellent controlled conditions for testing 
and debugging atmospheric models. Many experimenters including Phillips 
(1959), Gates (1962), Heburn (1972), Holloway, Spelman,and Manabe (1973), 
Elias (1973), and Mihok (1974) have used the Haurwitz stream function to 
generate initial conditions. 

It should be noted that the equations used in the model are not those 
of a non-divergent atmosphere. Phillips (1959) has noted that the 
presence of divergence in the barotropic atmosphere will slow the rate 
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of progression of the flow pattern especially for small values of the 
wave number m . 



The Haurwitz stream function is used by Phillips (1959) as the 
forcing function to obtain the geopotential from the non-linear balance 
equation 



V 2 <J> 



fV 2 4> - Vip • Vf - 2(|^) 2 + 2(- 



"sx 2 dy ! 



) 



( 21 ) 



Equation (21) 

_i_ ^ i d 2 <p 

a 2 cos 2 0 3A 2 



in spherical coordinates becomes 

1 3 , 0 34 >., 

+ cos 6 36 (cOS 6 36 )] 



-t [ — i — 

a 2 cos 2 0 



3 2 i/j 
3A 2 



1 3 / _ 0 3^m . 1 ^ 3f 

cos 0 30 (C ° S 9 30 )] + 2 30 30 

a 



- 2 ( 



2 a 
a cos 0 



^t ) 2 + 



3A30 J 



4 2 fl 

a cos 0 



d 2 lp d 2 lp 

3A 30 



( 22 ) 



The solution for 4> of the non-linear balance equation may be stated as 
<J> = a 2 A(0) + a 2 B(0) sin mA + a 2 C(0) (2 sin 2 mA - 1) (23) 



where <p is the geopotential perturbation 



* 

A(0) = |(2fi+ B) cos 2 0 + \(r^) 2 cos 2 ® 0[(m+l)cos 2 0 +(2m 2 -m-2)- 

a 



2m 



2 a 
cos 0 



] , 



( 24 ) 
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2(ft*B) ^2 



cos™ 0 [ (m^ + 2m+2) - (nri-1)^ cos^ 0] 




(25) 



and 



c(0) 




0 [(m+-l) cos^ 0 - (m+2)] 



(26) 



a 



This solution was obtained by Phillips (1959), Equation (23) provides 
a geopotential perturbation field. The height at a particular pressure 
level is obtained by adding the height perturbation from (23) to the 
mean height for that level obtained from the NACA standard atmosphere 
(Haltiner and Martin, 1957). 

The linear balance equation may be obtained by dropping the last 
two terms in (21) for Cartesian coordinates and (22) for spherical co- 
ordinates. The linear solution is given by 



(p =* F(6) + G(0) sin (mA - Vt) 



(27) 



where 




(28) 



and 




+ 2m~ + 2 
(nrt-1) (nrt-2) 



] sin(mX-vt) cos™ 0 




sin (mX - Vt) cos”*^ 0 



(29) 
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The solution to the linear balance equation produces initial fields 
which are in geostrophic equilibrium while the non-linear solution pro- 
duces fields which are in approximate gradient equilibrium and therefore 
more realistic. This is especially true near the equator where the 
geostrophic approximation is invalid. For this reason non-linear 
balancing is used in all experiments in this thesis except Experiment 3 
which involves a flow pattern of wave number 12. 

In Experiment 3, non-linear balancing results in aliasing, (Haltiner, 
1971), which distorts the initial fields. This problem was eliminated 
by the use of linear balancing. The effect of eliminating the non- 
linear terms is minimal because the wave amplitude used in the experi- 
ments is small. 



B. ANALYTIC WINDS 

Analytic zonal and meridional wind components are obtained from the 
derivatives 



and 



_ 1 jty 
a 90 



(30) 



1 dip 

a cos 0 9A 



(31) 



The results are 

u = - “[A*sin(mA-vt)cos m+ ^0 -mA sin(mA-Vt)cos m ^0sin^0-Ba^cos0] 
a 

(32) 

v = ~[A m sin0 cos m ^ 0 cos(mA - Vt) ] • (33) 

a 
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C. ANALYTIC TEMPERATURES 



Analytic temperatures, constant on each pressure surface, are 
generated in accordance with the National Advisory Committee for 
Aeronautics (NACA) standard atmosphere (Haltiner and Martin, 1957). 
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IV. RESULTS 



The experiments conducted for this thesis fall into two categories: 
(1) experiments comparing second order differencing versus mixed 
second and fourth order differencing as discussed in Section II-D-1 
and (2) experiments attempting to improve the forecast effectiveness 
of the model when there is flow over the poles. 

In all experiments presented in this thesis analytically derived 

fields are used as described in Section III. The variation in the 

character of the initial harmonic waves was accomplished by varying 

the specification of the wave number, phase speed and the amplitude 

* 

constant A which is part of the coefficient A given in Equation (19). 

The use of analytic fields resulted in the elimination of spurious 
inertial-gravity waves resulting from initial inconsistencies between 
the and <p fields. Thus, results directly related to the finite 

differencing of the primitive equations can be examined under con- 
trolled conditions without distortion due to initialization. 

Surface pressure analysis and forecasts are used to illustrate the 
output of the model. These analyses, for the Northern Hemisphere only, 
were obtained by interpolating the spherical grid onto the standard 
63 x 63 polar stereographic grid used by FNWC. The interpolation of 
the spherical data field for a point on the 63 x 63 grid is accomplished 
by using the Ayres 1 central difference formula, which is continuous in 
the first derivative. If the point lies within the border zone, a 
double linear interpolation is performed. Very little detail is lost, 
except at the pole as shall be shown in Experiment VIII. 
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Harmonic analysis (Jeffreys and Jeffreys, 1956) is performed on the 
surface pressure analysis and on the forecasts at three-hour intervals. 
The harmonic analysis reveals the wave numbers present ,around a latitude 
circle and their amplitudes and phases. 

A. SECOND ORDER DIFFERENCING VERSUS MIXED SECOND AND 

FOURTH ORDER DIFFERENCING 

Experiment I . In this experiment, wave number 6, a phase speed of 
+ 12° longitude per day and A =* 1000 are used. Chart A shows the initial 
surface pressure field. Chart B is the 48-hour forecast using second 
order differencing and Chart C shows the 48-hour forecast using the mixed 
differencing scheme. Note that the forecast fields retain their defini- 
tion and are not distorted by undesirable inertial gravity waves. 

Figure 3 shows the phase angles in degrees longitude as a function of 
latitude at 12-hour intervals to 48 hours for second order differencing. 
Figure 4 shows the same relationship for the mixed differencing scheme. 
Comparison of Figures 3 and 4 shows a considerable improvement in the 
forecast phase speed of the waves when the mixed differencing scheme is 
used. This is especially true as the latitude increases. Above 60° 
latitude the mixed scheme forecasts of phase speed are too fast while 
the second order forecasts indicate a greatly retarded phase speed. 

Experiment II . Wave number 9, a phase speed of + 12° longitude per 
day and A =* 0.1 are used in this experiment. Chart D shows the initial 
surface pressure field. Chart E is the 48-hour forecast using second 
order differencing. Chart F is the 48-hour forecast using the mixed 
differencing scheme. The forecast and initial fields are not distorted 
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by undesirable gravity waves. Figure 5 shows the phase angles in degrees 
longitude as a function of latitude at 12-hour intervals out to 48 hours 
for second order differencing. Figure 6 shows the same relationship for 
the mixed differencing scheme. Comparison of Figures 5 and 6 indicates 
a closer approximation to the true phase speed when the mixed scheme is 
used. Above 60° latitude the mixed scheme phase speeds are too fast while 
considerable slowing is evident for second order differencing. 

Experiment III . Wave number 12, a phase speed of + 12° longitude per 
day and A =* 5 x 10 ^ are used in this experiment. The initial surface 
pressure field is shown in Chart G while Charts H and I are the 48-hour 
forecast fields for the second order differencing and mixed differencing 
schemes, respectively. Note that the initial and forecast fields are 
not distorted by gravity waves. Figure 7 shows the phase angles in 
degrees longitude as a function of latitude at 12-hour intervals out to 
48 hours for second order differencing while Figure 8 shows the same 
relationship for the mixed differencing scheme. When the second order 
differencing is used there is very little wave propagation. The phase 
progression is very distorted and marked retrogression is noted above 
40° latitude. Figure 8 indicates a distinct improvement in the propa- 
gation of the waves when the mixed differencing scheme is employed. In 
Figure 8 the phase speed of the waves increases with latitude. However, 
at all latitudes the phase speeds are slower than for the nondivergent 
analytic solution. 

It should be noted that the solution to the linear balance equation 
was used to generate the initial height fields in this case. The use 
of the non-linear solution resulted in undesirable aliasing. 
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B. CROSS-POLAR FLOW 



McCollough (1974) remarked that "noise” developed in the polar 
regions as the forecast time progressed and that longer integrations 
indicated a need for a review of the polar region finite differencing 
and the method of handling the pole point. In order to examine the 
effectiveness of the model in forecasting for the polar regions, flow 
over the poles was generated analytically using wave number 1 and a 
phase speed of -94.4° longitude per day. After some experimentation 
it was decided to make A = 3 x 10^ yielding a maximum wind of 20 kts 
in the polar region. This method of generating cross-polar flow was 
taken from Holloway, Spelman and Manabe (1973). 

Chart J is the initial surface pressure analysis used for the 
flow over the pole experiments. Charts K-N illustrate surface pressure 
forecasts at 12-hour intervals out to 48 hours. All cross polar flow 
experiments use second order differencing. 

Examination of Chart J reveals that the isobars are evenly spaced 
and curve smoothly over the pole. However, the pattern of the isobars 
becomes increasingly distorted near the pole with progressive fore- 
casts as can be seen in Charts K-N. The objective of the subsequent 
experiments is to test various modifications to the model in order to 
improve the accuracy of forecasts in the polar region when cross-polar 
flow is present. 

Experiment IV . Based on the assumption that the distorted flow 
pattern in the polar regions was the result of errors in the derivatives 
near the pole, smoothing of the derivatives was increased by squaring 
the stability coefficient, S , defined in Equation (13). This 
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modification did not result in smoother forecasts in the polar regions 
and in fact the disturbance was slightly augmented as can be seen by a 
comparison of Chart L and Chart 0. 

Experiment V , In this experiment removal of the irregularities in 
the polar region was attempted by applying the Arakawa smoothing tech- 
nique directly to the main variables at three-hour intervals. The 
Arakawa smoothing was applied to tt, u, v, T and q on the ninth and 
tenth time steps after the Matsuno step. Charts P and Q are the 24-hour 
and 48-hour forecasts made with this technique. Comparison of the above 
charts with Charts L and N reveal a reduction in the amplitude of the 
disturbance at the pole and a more uniform flow at the higher latitudes. 

Experiment VI . The promising results of Experiment V made it seem 
probable that application of Arakawa smoothing to the main variables 
at every time step would totally remove the distortion of the isobars 
near the pole. Although this method gave somewhat improved results over 
the case where the main variables were not smoothed, the effect was 
still not as good as that achieved in Experiment V. 

Experiment VII . Holloway, Spelman and Manabe (1973) transformed 
vector components into a polar stereographic mapping prior to Fourier 
filtering in order to prevent excessive smoothing of flow over the pole. 
Their Fourier filtering method is similar to the Arakawa smoothing tech- 
nique except that the shorter waves are completely removed. In this 
experiment the main variables were smoothed every three hours, as in 
Experiment V; but prior to smoothing, the wind components u and v were 
transformed to a polar stereographic projection plane. The method of 
transformation is the same as used by Holloway, Spelman and Manabe. 



38 



This experiment resulted in the model becoming unstable after 27 
hours. 

Experiment VIII , At this point in the experimentation it was noted 
that the interpolation technique used to transform the surface pressure 
field from spherical coordinates to a 63 x 63 polar stereographic grid 
was not adequately handling the pole point. Hence, the interpolation 
scheme was bypassed at the pole and the pole value on the spherical 
grid was inserted directly onto the 63 x 63 grid which is centered at 
the pole. Chart R and S show the result of this modification of the 
output interpolation scheme. The flow appears to be smoother indicating 
that a portion of the distortion in the polar regions was due to a 
faulty output interpolation scheme. 

Experiment IX . This experiment tests the effect of using only 
averages around the 85° latitude circle for values needed at the poles. 
The modified output interpolation scheme (Experiment VIII) is included 
in the model for this experiment. The resulting 24-hour forecast is 
shown in Chart T which shows a considerable improvement over Chart L 
(the original cross-polar flow 24-hour forecast). However, there is 
little difference between the 24-hour forecast using just the modified 
output interpolation scheme (Experiment VIII, Chart S) and Chart T which 
includes the additional feature of using averages for polar parameters. 

Experiment X . This experiment consisted of a combination of the 
modified output interpolation scheme (Experiment VIII) and the applica- 
tion of Arakawa smoothing to the main variables at three-hour intervals 
(Experiment V) . Chart U shows that Arakawa smoothing of the main 
variables at three-hour * intervals reduces the distortion at the pole by 
a considerable amount. The improvement was larger than in Experiment V. 
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A 48-hour forecast was not produced in this experiment due to time 
limitations. 
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FIGURE 1. LOCATION OF VARIABLES. 

i is the column counter and j is the row counter. 
Point (1,1) is at 85s, 10E. 
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FIGURE 2. VERTICAL LAYERING 

In the above figure, sigma (0) is the dimensionless vertical 
coordinate, u and v are the zonal and meridional wind 
components, respectively, q is the specific humidity, <J> is 
the geopotential-, m is the terrain pressure and w is the 
vertical velocity (-a). 
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Figure 3. Phase angle (degrees longitude) vs latitude for second order 
differencing, wave number 6, phase speed 12°/day, A = 1000. (Latitudes 
with zero wave amplitude are not included.). 
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LAT I TUDE 




Figure 4. Phase angle (degrees longitude) vs latitude for mixed second 
and fourth order differencing for wave number 6, phase speed 12°/day, 

A = 1000. (Latitudes with zero wave amplitude are not included.) 
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Figure 5. Phase angle (degrees longitude) vs latitude for second order 
differencing, wave number 9, phase speed 12°/day, A = 0.1. (Latitudes 
with zero wave amplitude are not included.) 
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Figure 6. Phase angle (degrees longitude) vs latitude for mixed second 
and fourth order differencing, wave number 9, phase speed 12°/day, 

A » 0.1, (Latitudes with zero wave amplitude are not included.) 
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Figure 7. Phase angle (degrees longitude) vs latitude for second order 
differencing, wave number 12, phase speed l-2°/day, A = 5 x 10 . 

(Latitude with zero wave amplitudes are not included.) 
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Figure 8. Phase angle (degrees longitude) vs latitude for mixed second 
and fourth order differencing, wave number 12, phase speed 12°/day, 

A = 5 x 10~6 , (Latitudes with zero wave amplitudes are not included.) 
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Chart A. Initial surface pressure analysis, wave number 6, phase 
speed 12°/day f A = 1000, 
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Chart B. 48-hour surface pressure forecast second order differencing, 
wave number 6, phase speed 12°/day, A = 1000. 




Chart C, 48-hour surface pressure forecast, mixed second and fourth 
order differencing, wave number 6, phase speed 12°/day, A =* 1000, 




Chart D. Initial surface pressure analysis, wave number 9, phase 
speed 12°/day, A =* 0.1. 




Chart E. 48-hour surface pressure forecast, second order differencing, 
wave number 9, phase speed 12°/day, A = 0.1. 
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Chart F. 48-hour surface pressure forecast, mixed second and fourth 
order differencing, wave number 9, phase speed 12°/day, A = 0.1. 
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Chart G. Initial surface pressure analysis, wave number 12, phase 
speed 12°/day, A = 5.0 x 10 - ^ . 
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Chart H. 48-hour surface pressure forecast, second order differencing, 
wave number 12, phase speed 12°/day, A = 5.0 x 10 ^ . 
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Chart I. 48-hour surface pressure forecast, mixed second and fourth 
order differencing, wave number 12, phase speed 12°/day, A = 5.0 x 10 
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Chart J. Initial surface pressure analysis, wave number 1, phase 
speed -94.4°/day, A = 3.0 x 10^ . 
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Chart K. 12-hour surface pressure forecast, wave number 1, phase 
speed -9A.A°/day, A =* 3.0 x 10^ • 
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Chart L. 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3.0 x 10^ . 
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Chart M. 36-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3,0 x 10^ • 
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Chart N. 48-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3.0 x 10^ . 
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Chart 0. 24-hour surface pressure forecast, wave number 1, phase 
speed -94. 4° /day, A - 3.0 x 10^; Arakawa smoothing coefficient squared. 
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Chart P. 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3.0 x 10^; Arakawa smoothing applied to main 
variables at 3-hour intervals. 
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Chart Q. 48-hour surface pressure forecast, wave number 1, phase 
speed -94. 4° /day, A = 3.0 x 10^, Arakawa smoothing applied to main 
variables at 3-hour intervals. 
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Chart R. 12-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A =* 3.0 x 10^, modified output interpolation 
scheme. 
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Chart S* 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3,0 x 10^, modified output interpolation scheme. 
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Chart T . 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A * 3,0 x 10 7 , modified output interpolation scheme, 
averages used for polar values. 
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Chart U. 24-hour surface pressure forecast, wave number 1, phase 
speed -94.4°/day, A = 3.0 x 10 , modified output interpolation 
scheme, Arakawa smoothing applied to main variables at 3-hour 
intervals . 
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V. CONCLUSIONS 



The experiments comparing second order differencing and the mixed 
second and fourth order differencing scheme indicate that in general 
the mixed differencing scheme provides a much closer approximation to 
the true phase speed. Phase speeds using second order differencing 
tended to decrease and therefore become more inaccurate with increas- 
ing latitude. The mixed differencing scheme produced phase speeds 
which tended to increase with latitude therefore becoming more accurate. 
As the wavelength of the harmonic waves became shorter the mixed 
differencing scheme gave increasingly better phase speeds in com- 
parison with the second order phase speeds. This is especially evident 
in Experiment III where wave number 12 was used. It should be noted 
however that both schemes when compared only with themselves on the 
basis of wave number became steadily less accurate as the wave number 
increased. 

The distortion in the cross-polar flow was in part due to the 
interpolation scheme which transformed the spherical grid to a 63 x 63 
polar stereographic grid. A great deal of the remaining "noise" 
at the pole can be removed by applying Arakawa smoothing at three- 
hour intervals to the main variables. 

The application of Arakawa smoothing on a polar stereographic plane 
projection for the main vector variables caused the model to become 
unstable. Squaring the Arakawa smoothing coefficient did not have a 
beneficial effect on the cross-polar flow. Application of the Arakawa 



70 



smoothing technique to the main variables at every time step smoothed the 
cross-polar flow somewhat but was not as effective as the three-hour 
application* The use of 85° latitude averages for the polar quantities 
led to a small improvement in the forecast. 
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APPENDIX A 



FINITE DIFFERENCE EQUATIONS 



The finite difference equations given here are basically the same 
as those given in Mihok (1974) . Additions have been made to include 
special handling of the equations near the pole and the discussion of 
the determination of the gradient of terrain pressure has been changed. 
Friction, diffusion, and heating and moisture source terms are not 
included. In the following equations At is the time increment, AX 
is a five degree longitudinal grid increment and A0 is a 5 degree 
latitudinal grid increment, and Acr =* 0^2 is the vertical increment 
between layers. Refer to Figure 1 for the staggered spherical grid 
arrangement and subscripting. Figure 2 illustrates the vertical grid. 
It should be noted that in the following difference equations the k 
subscript on w indicates even sigma levels (0.8, 0.6, 0.4, 0.2) and 
for all other quantities indicates odd sigma levels (0.9, 0.7, 0.5, 

0.3, 0.1). 

Longitudinal momentum equation 
1. Left hand side 



8utt _ 




( 1 ) 



2At 



where 




+ IT. 



L, j-1 ,k 



) 



(l.D 



In the above Equation (1.1) the subscript 



L » i (i even) 



L - i-1 (j odd) 
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2. Right hand side 
a. Term one 



1 r 9 (uutt ) 9(uvtt cos 0 ) , 

a cos 6 1 9A 96 J 



- (u. + U. ,, . . ) U 7 r . . -(u. + U. . . , ) U 7 T. , /0 . , 

1 r ijk l+l,;)^ 1+1/2, j ,k ljk 1-1, J ,k y x-l/2,j > k 

a cos 0. 2 • 2AA 

J 



(u...+ u. _ .) V7T . . cos0 -(u + u. . „ . )vir. . COS0. . 

, k i,.i+2,k i,j-l,k n+1 ijk i,j-l,k .i-l , 

2 • 2A0 J 



( 2 ) 



where tt is computed as in (1.1) and 



_ l r/ - 



u7T i+l/2,j,k " 4^ U7T) ijk +(uTr) i+l,j,k +< ' u7T) L,j+l,j +(u7T) L,j-l,k ] (2,1) 



UTT 



. - /9 . , is obtained from (2.1) by incremental subscripting. 
1 “J_/ z, j ,ic 



W i,j+l,k " 4 t(v7T) ijk +(v7T) i,j+2,k +(vTr) L,j+l,k +(v7T) M,j+l,k ] (2,2) 



vTT . . _ , is obtained from (2.2) by incremental subscripting. 

1 y J“1 y ^ 

The subscripts 

L = i (j odd), 

L = i+1 (j even) , 

M =* i-1 (j odd), and 
M =* i (j even) . 
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When computing the above expressions around the 85° and 80° latitude 
circles any terra requiring an undefined polar quantity is ignored. 
For example 



( u * . i . , + U . ) 

i+ 1 >j .> k 

2 u7T i+l/2,j,k 



is not included in the computation at 85°S since the , 

1+1/2, j,k 

contains the undefined quantity ( u tt) t . . 

b . Term two 

w. (u. . + u ) - w. T (n . . ! + u. . , n ) 

8a " xj 1 2Aa 1 

where the vertical velocity is defined as 



Aa r/ 9Ti\ 

w. = w . . . i - — [ (~~) . . 

xjk i,j,k-l 7T 8t xj 



cos0 






. v-.,, ,COS0 -(vtt). . COS0. 

1 ' x+l^^k x.lk x». 1 +l»k , 1+1 x,.i-l,k . 1-1 



2AX 



2A0 



)] 



(4) 



and the local chance in terrain pressure as 



(— ) = 

ot'ij - 



-x 



Aa 



a cos0 . 

1 J 



1 2AX 



(V7T) . COS0 -(V7T). . - .COS0. 

x.i+j.k , 1+1 x,.i-l,k ,1-1 - 



2A0 



(5) 



The TT in the above equations is computed as in (1.1). The w terms 
are computed in the same manner as the TT terms. 
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When computing w and the local change in terrain pressure at the 
north pole for the thermodynamic and moisture equations the following 
expressions are used, as discussed in Section II-D-2, 



W pole,k 



W pole,k-l 



Aq r /3jT\ , 2A q 

7T . ^St^pole aA036 
pole r 



36 

Z (vS) 

i=l 



i,35,k ] 



( 6 ) 



and 



(— ) 
t '3t ; 




2Aa 

aA036 



36 



2 (v1T) i,35,k 

i=l 



(7) 



In order to compute these quantities at the south pole change the sign 
of (vrr) and the subscript j from 35 to 1 , 
c. Term three 



n (uirv) . .. tan0 . 

ttuv tanQ ilk i 

a a 

with 7T as in (1.1). 

d. Term four 

TTvf = f . (ttv) . . . 

J ijk 

with if as in (1.1). 

e. Term five 



( 8 ) 



(9) 



1 r RT in + Tt iii = 
a cos 0 ^ 3 A 3A 



1 << ’ i3k 2A> 1 ^- i ' ,k) ) 



a cos 0^ ijk v 3A ij ij 



( 10 ) 



where TT is obtained from (1.1) and T is computed in the same manner. 
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Based on Equation (10) of Section II-A the longitudinal derivative 



of terrain height may be expressed 



A* = 
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RT2AX 



[A* p - 4* 0 ] 



( 11 ) 



In order to obtain [A^ “ A$^] we use the hypsometric equation 

to form 



- A* 



RT [£n tt. - £n TT. 1 ] 
i l-l 



( 12 ) 



The right hand side of (12) can be modified to 



- Ad) = RT[£n7T .-£n7T . + i.n7T.- i,n7T. .. ] 
O 11 1 i-I 



(13) 



without changing the equality. 

Bringing T inside the parenthesis and computing it in a special way 
on either side of the velocity point yields 



[A<f> - A<f> ] . 

1 Y p T 0 ijk 



R[T LEFT,k (£n1T ij £n1T ij )+ T RlGHT (£n7T ij £n7r i-l,j^ 



(14) 



T LEFT anc * ^RIGHT are s P ec ^- a ^- mean temperatures derived over the five 
degree grid distance from the computational velocity point to the 
adjacent mass point either on the right or on the left. The following 
method is used to compute these temperatures* The subscript S indi- 
cates side (either right or left). 

(1) If (JinTT.. - JlnTT ) <0 or k =* 1 , (i.e. O = 0.9) then 
-LJ ^ 



T = t + tjl 

S,k A S,k 2 



< T S.k + l - T S,k> <* nll n - * n V 

(<lno k+r Jno k ) 



(15) 
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Equation (15) takes care of the situation where the pressure on the 
sigma surface at the computational point is less than the sigma surface 
pressure at an adjacent mass point. The mean change of temperature 
between the layers is multiplied by the logarithmic slope of the sigma 
surface • This is used to modify the temperature of the lower 

layer. The end result is to interpolate the value of the geopotential 
at the mass point to the same pressure surface as the computational 
point before integration takes place. 

(2) If (JlnTT . . - &n7T ) >0 or k = 5, (a = 0.1), then 



S,k 



_ (*ni..-£n7r s ) + qna k -in 

S,k-1 2 1 Jln(a k - In a p J 



(16) 



Equation (16) comes into play when the pressure on the sigma surface 
at the computational point is less than the pressure on the sigma sur- 
face at an adjacent mass point. The mean change of temperature between 
the layers is multiplied by the logarithmic slope of the pressure sur- 
face — . This is used to modify the upper layer. The end result 
is to interpolate the value of the geopotential at the mass point to the 
same pressure surface as the computational point before integration takes 
place. 

3 TT 

This computational scheme for * 7 ^* is intended to reduce truncation 
error when differencing over steep terrain. Even though a flat earth 
is used in this thesis the terrain pressure spatial derivatives using 
this scheme were allowed to remain. 
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Latitudinal momentum equation 



1. Left hand side 



3vtt . < iv >m - 



at 



2At 



( 17 ) 



where 7T is computed as in (1.1), 
2. Right hand side 

a. Term one 



1 r 3 (ttw ) 3(ttw cos 0) 1 

a cos0 . 1 3A 30 J 

3 



1 / V ljk +V i+l > ,i ,k^ UTr i+l/2, ,j ,k ^ijk^i-l, j t k^ U ^i-l/2,,j ,k 

a cos0 . 2 • 2AX 

J 



(v...+V. ... , )vTT. ... .COS0.J. - (v . 4-V . . _ . )VTT. . .. , C.OS0 . .. 

+ i,lk I,.i+2,k I,i+l,k 1+1 ilk 1,1-2, k i-l , 

2 • 2A0 J 



(18) 



where the utt and vtt terms are computed as in (2.1) and (2.2). 

b. Term two 



if- 



3 (wv) = " ^ijk^i. i .k+l^ijk^ W i , j , k-1 ^ V ij k +V i , ,j , k-1 1 

8a ij 1 2A a 2A a J 



(19) 



__ 3tt 

The w and (t^-) terms are computed in the same manner as they are 
for (3), the vertical acceleration term in the latitudinal momentum 
equation. 

c. Term three 



TTuutan 0 = (u7Tu) ijk tan 8 j 



( 20 ) 



where TV is computed as in (1.1). 



78 



d. Term four 



Tmf =* f . (ttu) , .. 

J ijk 



( 21 ) 



where tt is obtained as in (1.1) 
e. Term five 



- i [RT w + 77 w 1 5 



- — rRT (— )* + if — ^ l ?-3 ] 

a 1 ijk^e^ij ij 2A6 1 



(22) 



where 



L =* i (j even) 
L = i-1 (j odd) 



TT is obtained from (1.1) and T is computed in a similar manner. 
3 tt * 

The computation of is carried out in the same manner as the com- 

dt) 

putation of (-gy) for Equation (10) with the exception that the (-^g*) 
in this case is computed along a meridian instead of a latitude circle. 
Refer to Equations (10), (11), (12), (13), (14), (15), and (16). 
Thermodynamic equation 



1. Left hand side 

r,, ( TTT )?'!’, 1 ~ ( TTT )^. 1 

9(ttT) _ ink i.i k 



3t 



2At 



(23) 



This term is computed directly at the poles, 
2. Right hand side 
a. Term one 



1 . 3 (ttuT) 9 (ttvT cos 6 ) , __ 

a cos 0 3A 30 “a cos 0. 






2AA 



(Tvl,) L,W.k cose Hl - ^L.i-l.k^Vl . 

2A0 J 



(24) 
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where 



L =* i (j odd) 
L = i+1 (j even) 



77 is computed as in (1.1) and T is computed in a similar manner. At 

85°S (Tvtt) is undefined and not included in the difference equa- 

tion. At 85°N (Tvtt), . , - . is undefined and not included. The 

L,j+l,k 

computation of (24) at the poles is carried out as discussed in Section 
II-D-2 (Polar Finite Differencing). At the north pole the expression 
used for term one, (24), is 

36 

R.H.S. term 1 = - ^ £ <™0 1>35> k < 25 > 

i^l 

In order to compute this term at the south pole change the sign of 
(Tvtt) and the subscript j to 1 . 
b. Term two 



, / w. .. (T. . . , -+T. .. )-w. . . - (T. +T. . . .) 

3_(wT) _ r i.i k i,i, k+1 1.1 k I,i,k-1 n k i ,i,k-l 1 

71 3a ~ Vj 1 2Aa J 



(26) 



where Equations (4) and (5) are used to compute w . This term is 
computed directly at the poles, 
c. Term three 



RT f 

cuL' 



p 

RT 

C 



+ + v( 36 } cos e]) J = 



36" 



ijJi _ J- J . K U + n l A 

° k [ 2 k \ St « 



l . (u,l) n-i.i.k (ulT) iik _ lr ii (u i+i,i.k~ u iik ) 

' , a i 2AA 



2AA 



a cos0. 

+ cose. - cose. j,. . (v M +i , .k~ %J - 1 ,k) ]) 

OAQ J 1J 2A6 J y 



2A0 
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( 27 ) 



where 

L = i (j odd) 

L =» i+1 (j even) . 

3 TT 

w is computed by (A), (g^-) by (5), and tt by (1.1). Note a special 

scheme is used to compute the derivatives of terrain pressure in this 

equation in order to reduce truncation error. At 85°S - vtt t . - . and 

- v, . - . are undefined and not included in the difference equation. 
L,J-l,k 

At 85°N vtt . and v_ . , - . are undefined and not included. The 

L,J+1,J L,j+l,k 

computation of (27) at the poles is carried out as discussed in Section 
II-D-2 (Polar Finite Differencing). At the north pole the expression 
used for term three (27) is 



RT 

R.H.S. term 3 = P° le > - k . 

p \ 



W pole » k— 1^ ^pole 
2 






9tt 

at 



pole 



aA036 



^ (w pole,k + 

36 1 
5!1 [(vTr) i,35,k ” Tr i,35 V i,35,k ] jj 

i=*l 



(28) 



In order to compute this term at the south pole change the sign of 
(vtt) and (ttv) and change the subscript j to 1 . 

Moisture equation 



1. Left hand side 

, N n+1 f N n-1 
aqir _ (qTr) iik ~ (qTr) i1k 
at ' “ 2At 

This term is computed directly at the poles. 

2. Right hand side 
a. Term one 



(29) 
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1 , 9(uTiq) 3(vTTq cos 6 . 1 ^ u7T ^i+l, j ,k ^ U7T ^ijk 

a cos 0 3X 30 ^ “ a cos0 . 2AX 

J 



, ( « CT) L.1 + l.k C ° Se i + l - ^L.i-l.k^Vl 

2A0 



(30) 



where 

L = i (j odd) 

L = i+1 (j even) 



IT is computed as in (1.1) and q is computed in a similar manner. 

At 85°S (qvrT) . - , is undefined and at 85°N (qvir) _ , is 

L,J+1,K 

undefined. Undefined terms are not included in the difference equations. 
The computation of (30) at the poles is as discussed in Section II-D-2 
(Polar Finite Differencing) . At the north pole the expression used for 
term one, (30) , is 



36 

R.H.S. tenn 1 5 -^§55 ^ <5 vi >i,35,k: 

i=l 



(31) 



In order to compute this term at the south pole reverse the sign of 
(qvTT) and change the subscript j to 1 . 



b. Term two 



3 (wq) _ r W iik^ q i,j ,k+l +q ijk^ ~ W i, j ,k-l^ q ijk + q i,j,k-l^ .. 
30 lj 1 2A a J 



(32) 



where Equations (4) and (5) are used to compute w • This term is com- 
puted directly at the poles. 

Hydrostatic equation 

The hydrostatic equation is integrated to form the hypsometric 
equation 

A<J> * -RT A Jin a (33) 
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This is converted to a difference equation which is a modified form of 
that used by Kesel and Winninghoff (1972). 






i > j > k+1 



. Rp A±1 T j j k ( &nq k +&nTT i1 ) +T j t j , k +l ana k+l +&nTT jj ) 
9 ijk "V n £n 0, + £na, , + 2£mr_. J 



] 



(34) 



k+l 



iJ 



For the lowest layer (34) is modified to form 

♦iji ’ - EI iji < 35 > 

Both (34) and (35) are computed directly at the poles. 
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APPENDIX B 



TERMS USING FOURTH ORDER FINITE DIFFERENCING 
IN MIXED SECOND AND FOURTH ORDER SCHEME 



Williams (1972) has concluded that fourth order differencing used 
only in the advection terms and the divergence terra of the continuity 
equation will improve the phase speed of meteorological waves while 
keeping the necessary reduction in the time step to a minimum* The 
following equations are the advection terms and the divergence term in 
fourth order finite difference form as developed by Mihok (1974) for 
use in the global model. 

The advection term for the latitudinal momentum equation is 

1 j-4 r ^ijk^i+l, j »k^ U7r j+l/2 !> j ,k ^ U ijk +U j-1, j ? k^ U7r j-l/2. j *k 

a cos 0 . l 3 1 2 • 2AX 

J 

(U...+U. . , o I )v7T. , COS0 . - - (u..+u. . 0 . )V7T. . - COS0 . . 

+ . x J k ijJtlA. J+ 1 u ±L i zLJl J." 1 .} 

2 • 2A0 1 



1 / U ijk +U i+2 > j ,k^ u71 i+l,j ,k ^ U ijk +U i-2^j ,k^i-l, j ,k 

“ 3 1 2 • 4AA 



2 • 4A0 11 



( 1 ) 



where the urr. , , utt . . . , , vtt . ... , and vtt. . . , are 

1+1/2, j,k 1-1/2, j,k i, J+l,k i,j-l,k 

given by (2.1) and (2.2) in Appendix A. Note that 
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( 1 . 1 ) 



u7T i+l,j,k = 

"4 f (u 71 } i j k + (u 71 } i+ 2 , j , k + (u 71 } i+ 1 , j + 2 , k + (u 71 } i+1 , j - 2 , k 1 

utt. . , is obtained from (1.1) by incremental subscripting. 
1-1 >J 9 K 



^1,1+ 2,k = 

T [( '™ ) ijt +(vi) i-l > j+2,k +<vi) i+l,j+2,k +<v;) 



i,j+4,k 



] 



( 1 . 2 ) 



and VTT . - is obtained from (1.2) by incremental subscripting. 

1-9 J ^ 9 K 

If one substitutes v everywhere for u in (1), the result is the 
longitudinal momentum equation. 

The advection term for the thermodynamic equation is given by 



- , (Tutt) . , , - (Tutt) . . - 

1 1+1. i,k y i,nk 



a cos 0. l 3 
J 



2AA 



(Tvtt) COS0 - (Tutt) v cos0 

, 1,,1+l.k 1+1 i,i-l,k i-l -i 

2A0 ; 



-f( 






4AA 



(T*tvw]*)i . +2 k cos 8 m - (i*[vi?l*) lj1 . 2ik cose 1 . 2 

4A0 



( 2 ) 



where 7T is computed as in (1.1) of Appendix A and T is computed in 
a similar manner. Note that 
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( 2 . 1 ) 



T = — fT + T + T + T 1 

1+1, j,k 4 1 L, j ,k L,j+l,k A M, j ,k A L,j-l,k J 

T i-l,j,k = 4 [T N,j,k + T N,j+l,k + T 0,j,k + T N,j,k ] 

T i,j+2,k = 4 [T L,j+2,k + T L,j+3,k + T M,j+2,k + T L,j+l,k ] . 



( 2 . 2 ) 



(2.3) 



and T. . 9 , is obtained from (2.3) by incremental subscripting. 

In the above equations 

L =* i+1 (j odd) 

L = i+2 (j even) 

M = i+2 (j odd) 

M = i+2 (j even) 

N = i-1 (j odd) 

N =* i (j even) 

0 = i (j odd) 

0 =* i-1 (j even) 

_ * _ 

The expressions for [uir] can be obtained by substituting uTT in 

for T in the above equations. 

Quantities with a bar, e.g. T, are applicable at velocity points 

— * 

while starred quantities, e.g. T , are applicable at mass points. The 
advection term for the moisture equation is the same as (2) with q 
substituted for T . 
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The divergence term of the continuity equation is 



i <">i.,m.k <: ° 39 n-r<">i.i-i,k e ° se i-i 1 

a cos0. l 3 x 2AA 2A6 • 5 

J 

l r rUTr] i+l t i,k~ [u7T] i-l,,j,k [vTr] i t .i+2,k C ° s6 i+2~ [vTr] i,i-2,k COs6 i-2 1 , 

“ 3 X 4AA 4A0 11 



(3) 



where ir is given by (1.1) in Appendix A and the [ujr] quantities 
are computed as in the thermodynamic equation. 

Mihok (1974) blended the mixed second and fourth order scheme with 
the normal second order scheme at 75° latitude. This was done to avoid 
differencing over the poles. 
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